Method of modeling multi-mode degradation process and predicting remaining useful life

ABSTRACT

Disclosed is a method of modeling a multi-mode degradation process and predicting a remaining useful life, which belongs to the technical field of health management. The method comprises the following steps: firstly collecting degradation data of equal-interval sampling; performing change point detection for the degradation data; performing clustering with a degradation rate as a feature for degradation segments segmented by change points; establishing a degradation model comprising mode switching, wherein the mode switching is described by one continuous-time Markov chain; estimating a Hurst exponent in a degradation process by a quadratic variation method; estimating a state transition probability matrix of the Markov chain and both a drift term coefficient and a diffusion term coefficient in each mode by a maximum likelihood method respectively; obtaining an obeying distribution of a drift term under an influence of state switching in a period of time in future based on a Monte Carlo algorithm; and obtaining a distribution of a remaining useful life with a given threshold. The distribution of the remaining useful life of a system or equipment comprising a plurality of degradation modes is predicted more accurately.

TECHNICAL FIELD

The present disclosure belongs to the technical field of healthmanagement, and in particular to a method of modeling a multi-modedegradation process and predicting a remaining useful life.

BACKGROUND

Prediction of a remaining useful life of industrial equipment canprovide effective information for preparing a maintenance strategy and aproduction decision for equipment, thereby reducing losses resultingfrom an equipment failure and ensuring security and reliability of asystem.

Degradation process modeling is a critical procedure for predicting aremaining useful life. To obtain an accurate prediction result of theremaining useful life, a degradation model should conform to an actualdegradation situation as possible. At present, in most existing methods,it is assumed that there are no large working condition change andmaintenance operation with equipment within its entire useful lifeperiod. However, in an actual industrial production, a plurality ofworking conditions may exist in an equipment running process, and theequipment will be subjected to regular or condition-based repair ormaintenance activities at the same time. These activities may affect adegradation process of the equipment, thereby presenting a plurality ofmodes in the degradation process. At present, a method of predicting aremaining useful life for automatically identifying a plurality ofdegradation modes still does not emerge.

Modeling of a multi-mode degradation process and prediction of aremaining useful life mainly have the following difficulties: firstly,it is required to identify the number of degradation modes and thedegradation model in each mode according to historical data sincedegradation mode switching does not have a label; secondly, it isdifficult to obtain an analyzed first hitting time distribution sincethe degradation model includes a fractal Brownian motion which isneither a Markov process nor a semi-martingale; thirdly, a futuredegradation mode switching condition is unknown, and thus it is requiredto consider possible degradation mode switching in future during theprediction of the remaining useful life.

SUMMARY

To solve the above technical problems in the prior art, the presentdisclosure provides a method of modeling a multi-mode degradationprocess and predicting a remaining useful life, which is reasonable indesign and overcomes the disadvantages of the prior art, producing agood effect.

To achieve the above objects, the following technical solution isadopted in the present disclosure.

A method of modeling a multi-mode degradation process and predicting aremaining useful life includes the following steps:

at step 1, collecting degradation data x₀, x₁, x₂, . . . , x_(k) ofequipment at equal-interval sampling moments t₀, t₁, t₂, . . . , t_(k)respectively, where a sampling interval is τ, and the number of samplingis k;

at step 2, detecting slope change points, denoted as γ₁, γ₂, . . . γ_(j)and γ_(j+1) . . . , of a historical degradation process according to achange point detection method;

at step 3, obtaining a degradation segment by taking the points γjD andγ_(j+1) obtained at step 2 as endpoints, calculating a slope η_(γ)_(j+1) of the degradation segment based on the following formula, andtaking the slope η_(γ) _(j+1) as a feature value of the j-th degradationsegment;

$\eta_{\gamma_{j + 1}} = \frac{\sum\limits_{j = i_{\gamma}}^{i_{\gamma + 1}}\left( {x_{j + 1} - x_{j}} \right)}{i_{\gamma + 1} - 1}$

calculating a local density ρ_(j) of the feature value of eachdegradation segment, and calculating a minimum distance δ_(j) of thefeature value greater than the local density, where the local densityρ_(j) is calculated based on the formula (1):

$\begin{matrix}{{\rho_{j} = {\sum\limits_{j}{\chi \left( {d_{ji} - d_{c}} \right)}}},} & (1)\end{matrix}$

where d_(c) is a truncation distance, d_(ji)=|η_(i)−η_(j)|, and afunction χ(·) is defined as follows:

$\begin{matrix}{{{\chi (a)} = \left\{ \begin{matrix}{1,{a < 0}} \\{0,{a \geq 0}}\end{matrix} \right.};} & (2)\end{matrix}$

and

calculating the minimum distance δ_(j) based on the following formula(3):

δ_(j)=min_(i:ρ) _(i) _(>ρ) _(j) d _(ji)  (3);

at step 4, if ρ_(j) and δ_(j) are greater than corresponding thresholdsrespectively, a line segment obtained with the slope change points γ_(j)and γ_(j+1) as endpoints being one clustering center; according to thismethod, denoting the number of the obtained clustering centers as N,that is, clustering the line segments segmented by the slope changepoints as N categories according to the slopes of the line segments,denoting a degradation mode of a sampling moment u as Φ(u), lettingφ_(i)=Φ(t_(i)), and denoting time points of changes of the degradationmode as c₁, c₂, . . . ;

at step 5, establishing a degradation model based on the formula (4):

X(t)=X(0)+∫₀ ^(t)λ[Φ(u)]du+σ _(H) B _(H)(t)  (4),

where X(0) refers to an initial value of a degradation process, λ[Φ(u)]refers to a drift term coefficient; if λ(i)˜N(μ_(λ) _(i) , σ_(λ) _(i)²), σ_(H) refers to a diffusion term coefficient, B_(H) (t) refers to astandard fractal Brownian motion, and the degradation mode Φ(u) is onecontinuous-time Markov chain with a transition probability matrix beingQ;

${Q = \begin{bmatrix}{- q_{1}} & q_{12} & \ldots & q_{1N} \\q_{21} & {- q_{2}} & \ldots & q_{2N} \\\vdots & \vdots & \ddots & \vdots \\q_{N1} & q_{N2} & \ldots & {- q_{N}}\end{bmatrix}};$

at step 6, estimating q_(j) and q_(ij) in the transition probabilitymatrix Q of the continuous-time Markov chain according to φ₀, φ₁, φ₂, .. . φ_(k) based on the formulas (5) and (6):

$\begin{matrix}{{q_{j} = \frac{m_{j}}{\sum\limits_{i = 1}^{m,}v_{i}^{(j)}}};} & (5) \\{q_{ij} = {\frac{m_{ij}}{m_{i}}{q_{i}.}}} & (6)\end{matrix}$

where m_(j) refers to a number of times that the degradation mode hitsand stays in the j-th mode before the moment t_(k), m_(ij) refers to anumber of times that the degradation mode transits from the i-th mode tothe j-th mode before the moment t_(k), and v_(i) ^((j)) refers to a staytime of the degradation mode hitting the j-th mode at the i-th time;

at step 7, estimating a Hurst exponent H of the degradation processbased on the formula (7):

$\begin{matrix}{H = {\frac{1}{2}\log_{2}\frac{E\left\lbrack \left( {\sum\limits_{j = 1}^{p}{\theta_{j}x_{2{({i + j})}}}} \right)^{2} \right\rbrack}{E\left\lbrack \left( {\sum\limits_{j = 1}^{p}{\theta_{j}x_{i + j}}} \right)^{2} \right\rbrack}}} & (7)\end{matrix}$

where θ₁, θ₂, . . . , θ_(p) refer to wavelet decomposition high-passfilter coefficients of Symlets wavelet function, p refers to a number oforder of a vanishing moment of the wavelet function, and E(·) refers toa mathematic expectation;

at step 8, estimating an estimation value λ_(i) of the drift termcoefficient λ[Φ(u)] in the degradation mode of the i-th segment based onthe formula (8) respectively:

$\begin{matrix}{{\lambda_{i} = \frac{{\overset{\sim}{x}}_{c_{i}:c_{i + 1}}^{T}Q_{{\overset{\sim}{x}}_{c_{i}:c_{i + 1}}}^{- 1}I_{i}}{\tau \; I_{i}^{T}Q_{{\overset{\sim}{x}}_{c_{i}:c_{i + 1}}}^{- 1}I_{i}}},} & (8)\end{matrix}$

where I_(i) refers to one c_(i+1)−c_(i)-dimension column vector, eachelement thereof is 1, {tilde over (x)}_(c) _(i) _(:c) _(i+1) =[x_(c)_(i) ₊₁−x_(c) _(i) , x_(c) _(i) ₊₂−x_(c) _(i) ₊₁, . . . , x_(c) _(i+1)−x_(c) _(i+1) ₊₁]^(T),

$Q_{{\overset{\sim}{x}}_{c_{i}:c_{i + 1}}}$

refers to one c_(i+1)−c_(i)-dimension covariance matrix, and the elementin the i-th row and the j-th column thereof is½[|i−j+1|^(2H)τ^(2H)+|i−j−1|^(2H)τ^(2H)−2|i−j|^(2H)τ^(2H)];

at step 9, estimating an expectation μ_(λ) _(j) and a variance σ_(λ)_(j) of the drift term coefficient λ[Φ(u)] in each degradation modebased on the formulas (9) and (10) respectively:

$\begin{matrix}{{\mu_{\lambda_{j}} = \frac{\sum\limits_{i = 1}^{m_{j}}\lambda_{j,i}}{m_{j}}};} & (9) \\{{{\sigma_{\lambda}}_{j} = \sqrt{\frac{\sum\limits_{i = 1}^{m_{j}}\left( {\lambda_{j,i} - \mu_{\lambda_{j}}} \right)^{2}}{m_{j}}}},} & (10)\end{matrix}$

where λ_(j,i) refers to an estimation value of the drift termcoefficient obtained when the j-th degradation mode is hit at the i-thtime;

at step 10, estimating a diffusion coefficient σ_(H) based on theformula (11):

$\begin{matrix}{{\sigma_{H} = \sqrt{\frac{\left( {{\overset{\sim}{x}}_{0:k}^{T} - {\Lambda\tau}} \right)^{T}{Q_{{\overset{\sim}{x}}_{0:k}}^{- 1}\left( {{\overset{\sim}{x}}_{0:k}^{T} - {\Lambda\tau}} \right)}}{k}}},} & (11)\end{matrix}$

where Λ^(T)=[λ₁I₁ ^(T), λ₂I₂ ^(T), . . . , λ_(m)I_(m) ^(T)], {tilde over(x)}_(0:k)=[x₁−x₀, x₂−x₁, . . . , x_(k)−x_(k- 1)]^(T),Q_({tilde over (x)}) _(0:k) refers to one k-dimension covariance matrix,and the element in the i-th row and the j-th column thereof is½[|i−j+1|^(2H)τ^(2H)+|i−j−1|^(2H)τ^(2H)−2|i−j|^(2H)τ^(2H)];

at step 11, letting Ω[Φ(t_(k)),l_(k)]=∫_(t) _(k) ^(t) ^(k) ^(+l) ^(k)λ[Φ(u)]du, and obtaining a numerical distribution ƒ_(Ω[Φ(t) _(k) _(),l)_(k) _(]) of Ω[Φ(t_(k)),l_(k)] by Monte Carlo method; and

at step 12, for a given failure threshold ω, an approximate distributionof a first hitting time of the degradation process being:

$\begin{matrix}{{{f_{l}\left( l_{k} \right)} = {\sum\limits_{i = 1}^{N}{\int_{\lambda_{\min}l_{k}}^{\lambda_{\max}l_{k}}{{p_{{\Phi {(t_{k})}}i}\left( l_{k} \right)}{f_{\Omega {({{\Phi {(t_{k})}},l_{k}})}}(s)}\frac{\sigma \left( l_{k} \right)}{\sqrt{2\pi {\sigma (0)}{\int_{0}^{l_{k}}{{\sigma (t)}dt}}}}\left\{ {\frac{\omega - x_{k} - s}{\int_{0}^{l_{k}}{{\sigma (t)}dt}} + \frac{\lambda (i)}{\sigma \left( l_{k} \right)}} \right\} \times \exp \left\{ {- \frac{\left( {\omega - x_{k} - s} \right)^{2}}{2{\sigma (0)}{\int_{0}^{l_{k}}{{\sigma (t)}dt}}}} \right\} {ds}}}}},} & (12)\end{matrix}$

where λ_(min)=min{λ(1), λ(2), . . . , λ(N)}, λ_(max)=max{λ(1), λ(2), . .. , λ(N)}, p_(Φ(t) _(k) _()i)(l_(k))=P{{|Φ(t_(k)+l_(k))}=i|Φ(t_(k))},and

$\begin{matrix}{{{\sigma (t)} = {\sigma_{H}\begin{Bmatrix}{{\overset{\lfloor\frac{t}{\tau}\rfloor}{\underset{i = 1}{\Sigma}}\left\lbrack {\int_{{({i - 1})} \cdot \tau}^{i \cdot \tau}{c_{H}s^{\frac{1}{2}}{\int_{{\lfloor\frac{t}{\tau}\rfloor} \cdot \tau}^{{\lfloor\frac{t + \tau}{\tau}\rfloor} \cdot \tau}{\left( {u - s} \right)^{H - \frac{3}{2}}u^{H - \frac{1}{2}}{duds}}}}} \right\rbrack}^{2} +} \\\left\lbrack {\int_{{\lfloor\frac{t}{\tau}\rfloor} \cdot \tau}^{{\lfloor\frac{t + \tau}{\tau}\rfloor} \cdot \tau}{c_{H}s^{\frac{1}{2}}{\int_{s}^{{\lfloor\frac{t + \tau}{\tau}\rfloor} \cdot \tau}{\left( {u - s} \right)^{H - \frac{3}{2}}u^{H\frac{1}{2}}{duds}}}}} \right\rbrack^{2}\end{Bmatrix}^{\frac{1}{2}}}}\mspace{20mu} {{c_{H} = \sqrt{\frac{2H{\Gamma \left( {\frac{3}{2} - H} \right)}}{{\Gamma \left( {\frac{1}{2} + H} \right)}{\Gamma \left( {2 - {2H}} \right)}}}},}} & {(13);}\end{matrix}$

where Γ(·) refers to Gamma function;

Preferably, step 2 specifically including the following steps:

at step 2.1, initializing a parameter, letting γ₁=0, i=1, and i_(γ)=1,and selecting a minimum interval mτ between two change points and athreshold ω_(β) of change point detection;

at step 2.2, calculating a slope of a degradation segment from aprevious change point to the i-th point; if i−i_(γ)>m, calculating theslope η_(i) of the current segment based on the formula (14), andletting i=i+1:

$\begin{matrix}{{\eta_{i} = \frac{\overset{i - 1}{\sum\limits_{j = i_{\gamma}}}\left( {x_{j + 1} - x_{j}} \right)}{i - i_{\gamma} - 1}};} & (14)\end{matrix}$

at step 2.3, calculating a change point detection index β(i) of the i-thpoint based on the formula (15):

$\begin{matrix}{{{\beta (i)} = \frac{{\sum\limits_{j = 1}^{m}\left( {x_{i} + {\eta_{i}j\; \tau} - x_{j}} \right)}}{\sqrt{1 + \eta^{2}}}};} & (15)\end{matrix}$

at step 2.4, determining whether the change point detection index β(i)of the i-th point exceeds the threshold ω_(β); if β(i)>ω_(β), x_(i)being a change point, and letting i_(γ)=i_(γ)+1 and γ_(i) _(γ) =i; and

at step 2.5, if i≤k, letting i=i+1, and then, performing step 2.2.

Preferably, step 11 specifically includes the following steps:

at step 11.1, selecting the number n of Monte Carlo samples, andinitializing parameters i=1 and v_(i,j)=φ_(k);

at step 11.2, generating n random numbers r_(j) obeying a uniformdistribution on [0,1];

at step 11.3, for the j-th Monte Carlo sample sequence, letting i=i+1and v_(i+1,j)=s, where s satisfies

${{\sum\limits_{w = 1}^{s - 1}{p_{v_{i,j},w}(\tau)}} < r_{f} \leq {\sum\limits_{w = 1}^{s}{p_{v_{i,j},w}(\tau)}}},$

and p_(v) _(i,j) _(,w)(τ) refers to a probability that the degradationmode transforms from the v_(i,j)-th mode into the w-th mode over time τ;if i<l_(k)/τ, returning to step 11.2; otherwise, performing step 11.4;

at step 11.4, calculating a total time length that each Monte Carlosample sequence stays in each mode within a time interval (t_(k),t_(k)+l_(k)), and denoting the length as

${S_{s,j} = {\sum\limits_{i = 1}^{l_{k}/\tau}{I\left( {v_{i,j} = s} \right)}}};$

and

at step 11.5, calculating a numerical distribution ƒ_(Ω(Φ(t) _(k) _(),l)_(k) ₎(x) of Ω[Φ(t_(k)),l_(k)] based on the formula (16):

$\begin{matrix}{{{f_{\Omega {({{\Phi {(t_{k})}},l_{k}})}}(x)} = {\frac{\underset{j = 1}{\sum\limits^{n}}{I\left( {{\overset{N}{\sum\limits_{s = 1}}{S_{s,j}{\lambda (s)}}} = x} \right)}}{n}\mspace{20mu} {when}}}{{{\sum\limits_{s = 1}^{N}{S_{s,j}{\lambda (s)}}} = x};}} & (16)\end{matrix}$

is established,

${{I\left( {{\sum\limits_{s = 1}^{N}{S_{s,j}{\lambda (s)}}} = x} \right)} = 1};$

otherwise,

${{I\left( {{\sum\limits_{s = 1}^{N}{S_{s,j}{\lambda (s)}}} = x} \right)} = 0}.$

Beneficial effects brought by the present disclosure are describedbelow.

A method of identifying a degradation mode from historical degradationdata in the present disclosure is applicable to a system with stateswitching, an environment change or a maintenance operation, and thus iscloser to an actual situation than the prior art. Different degradationmodes are obtained by identification according to the historical data, adegradation model that contains state switching and is based on afractal Brownian motion is established, and the switching of thedegradation mode is described through one continuous-time Markov chain.A state transition probability of the Markov chain and both a driftcoefficient and a diffusion coefficient in the degradation model areestimated by a maximum likelihood method respectively. In the presentdisclosure, to obtain the first hitting time of the degradation process,a numerical distribution of a diffusion term in a future time segment isfirstly obtained by a Monte Carlo method, and then, the first hittingtime of the degradation process based on the fractal Brownian motion isconverted into a time when a standard Brownian motion firstly hits atime-varying threshold containing an uncertainty through one time-spacetransformation, thereby obtaining a distribution of the remaining usefullife of the degradation process.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flowchart illustrating a method of modeling a degradationprocess and predicting a remaining useful life according to an exampleof the present disclosure.

FIG. 2 is a flowchart illustrating a change point detection methodaccording to an example of the present disclosure.

FIG. 3 is a flowchart illustrating a Monte Carlo method according to anexample of the present disclosure.

FIG. 4 is schematic diagram illustrating a temperature degradation curveof a furnace wall of a blast furnace according to an example of thepresent disclosure.

FIG. 5 is schematic diagram illustrating change point detection resultsof a temperature degradation curve of a furnace wall of a blast furnacein the first 1200 days according to an example of the presentdisclosure.

FIG. 6 is a schematic diagram illustrating degradation modeidentification results of a temperature degradation curve of a furnacewall of a blast furnace in the first 1200 days according to an exampleof the present disclosure.

FIG. 7 is a schematic diagram illustrating a prediction result of aremaining useful life of a furnace wall of a blast furnace according toan example of the present disclosure.

DETAILED DESCRIPTION OF THE EMBODIMENTS

The present disclosure will be further described below in detail incombination with accompanying drawings and specific examples.

A method of modeling a multi-mode degradation process and predicting aremaining useful life includes the following steps with a flow as shownin FIG. 1:

at step 1, collecting degradation data x₀, x₁, x₂, . . . , x_(k) ofequipment at equal-interval sampling moments t₀, t₁, t₂, . . . , t_(k)respectively, where a sampling interval is τ, and the number of samplingis k;

at step 2, detecting slope change points, denoted as γ₁, γ₂, . . . , ofa historical degradation process according to a change point detectionmethod (with a flow as shown in FIG. 2);

at step 3, obtaining a degradation segment by taking the points γ_(j)and γ_(j+1) obtained at step 2 as endpoints, calculating a slope η_(γ)_(j+1) , of the degradation segment based on the following formula, andtaking the slope η_(γ) _(j+1) as a feature value of the j-th degradationsegment;

$\eta_{\gamma_{j + 1}} = \frac{\underset{j = i_{\gamma}}{\sum\limits^{i_{\gamma + 1}}}\left( {x_{j + 1} - x_{j}} \right)}{i_{\gamma + 1} - 1}$

calculating a local density ρ_(j) of a feature value of each degradationsegment, and calculating a minimum distance δ_(j) of a feature valuegreater than the local density, where the local density ρ_(j) iscalculated based on the formula (1):

$\begin{matrix}{{\rho_{j} = {\sum\limits_{i}{\chi \left( {d_{ji} - d_{c}} \right)}}},} & (1)\end{matrix}$

where d_(c) is a truncation distance, d_(ji)=|η_(i)−η_(j)|, and afunction χ(·) is defined as follows:

$\begin{matrix}{{{\chi (a)} =}\left\{ {\begin{matrix}{1,{a < 0}} \\{0,{a \geq 0}}\end{matrix};} \right.} & (2)\end{matrix}$

and

calculating the minimum distance δ_(j) based on the following formula(3):

δ_(j)=min_(i:ρ) _(i) _(>ρ) _(j) d _(ji)  (3);

at step 4, if ρ_(j) and δ_(j) are greater than corresponding thresholdsrespectively, a line segment obtained with the slope change points γ_(j)and γ_(j+1) as endpoints being one clustering center; according to thismethod, denoting the number of the obtained clustering centers as N,that is, clustering the line segments segmented by the slope changepoints as N categories according to the slopes of the line segments,denoting a degradation mode of a sampling moment u as Φ(u), lettingφ_(i)=Φ(t_(i)), and denoting time points of changes of the degradationmode as c₁, c₂, . . . ;

at step 5, establishing a degradation model based on the formula (4):

X(t)=X(0)+∫₀ ^(t)λ[Φ(u)]du+σ _(H) B _(H)(t)  (4),

where X(0) refers to an initial value of a degradation process, λ[Φ(u)]refers to a drift term coefficient; if λ(i)˜N(μ_(λ) _(i) , σ_(λ) _(i)²), σ_(H) refers to a diffusion term coefficient, B_(H) (t) refers to astandard fractal Brownian motion, and the degradation mode Φ(u) is onecontinuous-time Markov chain with a transition probability matrix beingQ;

${Q = \begin{bmatrix}{- q_{1}} & q_{12} & \ldots & q_{1N} \\q_{21} & {- q_{2}} & \ldots & q_{2N} \\\vdots & \vdots & \ddots & \vdots \\q_{N1} & q_{N2} & \ldots & {- q_{N}}\end{bmatrix}};$

at step 6, estimating q_(j) and q_(ij) in the transition probabilitymatrix Q of the continuous-time Markov chain according to φ₀, φ₁, φ₂, .. . , φ_(k) based on the formulas (5) and (6):

$\begin{matrix}{{q_{j} = \frac{m_{j}}{\sum\limits_{i = 1}^{m_{j}}v_{i}^{(j)}}};} & (5) \\{{q_{ij} = {\frac{m_{ij}}{m_{i}}q_{i}}},} & (6)\end{matrix}$

where m_(j) refers to a number of times that the degradation mode hitsand stays in the j-th mode before the moment t_(k), m_(ij) refers to anumber of times that the degradation mode transits from the i-th mode tothe j-th mode before the moment t_(k), and v_(i) ^((j)) refers to a staytime of the degradation mode hitting the j-th mode at the i-th time;

at step 7, estimating a Hurst exponent H of the degradation processbased on the formula (7):

$\begin{matrix}{{H = {\frac{1}{2}\log_{2}\frac{E\left\lbrack \left( {\sum\limits_{j = 1}^{p}{\theta_{j}x_{2{({i + j})}}}} \right)^{2} \right\rbrack}{E\left\lbrack \left( {\sum\limits_{j = 1}^{p}{\theta_{j}x_{i + j}}} \right)^{2} \right\rbrack}}},} & (7)\end{matrix}$

where θ₁, θ₂, . . . , θ_(p) refer to wavelet decomposition high-passfilter coefficients of Symlets wavelet function, p refers to a number oforder of a vanishing moment of the wavelet function, and E(·) refers toa mathematic expectation;

at step 8, estimating an estimation value λ_(i) of a drift termcoefficient λ[Φ(u)] in the degradation mode of the i-th segment based onthe formula (8) respectively:

$\begin{matrix}{{\lambda_{i} = \frac{{\overset{\sim}{x}}_{c_{i}:c_{i + 1}}^{T}Q_{{\overset{\sim}{x}}_{c_{i}:c_{i + 1}}}^{- 1}I_{i}}{\tau I_{i}^{T}Q_{{\overset{\sim}{x}}_{c_{i}:c_{i + 1}}}^{- 1}I_{i}}},} & (8)\end{matrix}$

where I_(i) refers to one c_(i+1)−c_(i)-dimension column vector, eachelement thereof is 1, {tilde over (x)}_(c) _(i) _(:c) _(i+1) =[x_(c)_(i) ₊₁−x_(c) _(i) , x_(c) _(i) ₊₂−x_(c) _(i) ₊₁, . . . , x_(c) _(i+1)−x_(c) _(i+1) ₊₁]^(T),

$Q_{{\overset{\sim}{x}}_{c_{i}:c_{i + 1}}}$

refers to one c_(i+1)−c_(i)-dimension covariance matrix, and the elementin the i-th row and the j-th column thereof is½[|i−j+1|^(2H)τ^(2H)+|i−j−1|^(2H)τ^(2H)−2|i−j|^(2H)τ^(2H)];

at step 9, estimating an expectation μ_(λ) _(j) and a variance σ_(λ)_(j) of the drift term coefficient λ[Φ(u)] in each degradation modebased on the formulas (9) and (10) respectively:

$\begin{matrix}{{\mu_{\lambda_{j}} = \frac{\sum\limits_{i = 1}^{m_{j}}\lambda_{j,i}}{m_{j}}};} & (9) \\{{\sigma_{\lambda_{j}} = \sqrt{\frac{\sum\limits_{i = 1}^{m_{j}}\left( {\lambda_{j,i} - \mu_{\lambda_{j}}} \right)^{2}}{m_{j}}}},} & (10)\end{matrix}$

where λ_(j,i) refers to an estimation value of the drift termcoefficient obtained when the j-th degradation mode is hit at the i-thtime;

at step 10, estimating a diffusion item coefficient σ_(H) based on theformula (11):

$\begin{matrix}{{\sigma_{H} = \sqrt{\frac{\left( {{\overset{\sim}{x}}_{0:k}^{T} - {\Lambda \tau}} \right)^{T}{Q_{{\overset{\sim}{x}}_{0:k}}^{- 1}\left( {{\overset{\sim}{x}}_{0:k}^{T} - {\Lambda \tau}} \right)}}{k}}},} & (11)\end{matrix}$

where Λ^(T)=[λ₁I₁ ^(T), λ₂I₂ ^(T), . . . , λ_(m)I_(m) ^(T)], {tilde over(x)}_(0:k)=[x₁−x₀, x₂−x₁, . . . , x_(k)−x_(k- 1)]^(T),Q_({tilde over (x)}) _(0:k) refers to one k-dimension covariance matrix,and the element in the i-th row and the j-th column thereof is½[|i−j+1|^(2H)τ^(2H)+|i−j−1|^(2H)τ^(2H)−2|i−j|^(2H)τ^(2H)];

at step 11, letting Ω[Φ(t_(k)),l_(k)]=∫_(t) _(k) ^(t) ^(k) ^(+l) ^(k)λ[Φ(u)]du, and obtaining a numerical distribution ƒ_(Ω[Φ(t) _(k) _(),l)_(k) _(]) of Ω[Φ(t_(k)),l_(k)] by Monte Carlo method; and

at step 12, for a given failure threshold ω, an approximate distributionof a first hitting time of the degradation process being:

$\begin{matrix}{{{f_{l}\left( l_{k} \right)} = {\sum\limits_{i = 1}^{N}{\int_{\lambda_{\min}l_{k}}^{\lambda_{\max}l_{k}}{{p_{{\Phi {(t_{k})}}i}\left( l_{k} \right)}{f_{\Omega {({{\Phi {(t_{k})}},l_{k}})}}(s)}\frac{\sigma \left( l_{k} \right)}{\sqrt{2\pi {\sigma (0)}{\int_{0}^{l_{k}}{{\sigma (t)}dt}}}}\left\{ {\frac{\omega - x_{k} - s}{\int_{0}^{l_{k}}{{\sigma (t)}dt}} + \frac{\lambda (i)}{\sigma \left( l_{k} \right)}} \right\} \times \exp \left\{ \frac{\left( {\omega - x_{k} - s} \right)^{2}}{2{\sigma (0)}{\int_{0}^{l_{k}}{{\sigma (t)}dt}}} \right\} {ds}}}}},} & (12)\end{matrix}$

where λ_(min)=min{λ(1), λ(2), . . . , λ(N)}, λ_(max)=max{λ(1), λ(2), . .. , λ(N)}, p_(Φ(t) _(k) _()i)(l_(k))=P{{|Φ(t_(k)+l_(k))}=i|Φ(t_(k))},and

$\begin{matrix}{{{\sigma (t)} = {\sigma_{H}\begin{Bmatrix}{{\sum\limits_{i = 1}^{\lfloor\frac{t}{\tau}\rfloor}\left\lbrack {\int_{{({i - 1})} \cdot \tau}^{i \cdot \tau}{c_{H}s^{\frac{1}{2}}{\int_{{\lfloor{\frac{t}{\tau}}\rfloor} \cdot \tau}^{{\lfloor\frac{t + \tau}{\tau}\rfloor} \cdot \tau}{\left( {u - s} \right)^{H - \frac{3}{2}}u^{H - \frac{1}{2}}{duds}}}}} \right\rbrack^{2}} +} \\\left\lbrack {\int_{{\lfloor\frac{t}{\tau}\rfloor} \cdot \tau}^{{\lfloor\frac{t + \tau}{\tau}\rfloor} \cdot \tau}{c_{H}s^{\frac{1}{2}}{\int_{s}^{{\lfloor\frac{t + \tau}{\tau}\rfloor} \cdot \tau}{\left( {u - s} \right)^{H - \frac{3}{2}}u^{H - \frac{1}{2}}{duds}}}}} \right\rbrack^{2}\end{Bmatrix}^{\frac{1}{2}}}};} & \left. 13 \right) \\{\mspace{79mu} {{c_{H} = \sqrt{\frac{2H{\Gamma \left( {\frac{3}{2} - H} \right)}}{{\Gamma \left( {\frac{1}{2} + H} \right)}{\Gamma \left( {2 - {2H}} \right)}}}},}} & \;\end{matrix}$

where Γ(·) refers to Gamma function;

Step 2 specifically includes the following steps (with the flow as shownin FIG. 2):

at step 2.1, initializing a parameter, letting γ₁=0, i=1, and i_(γ)=1,and selecting a minimum interval mτ between two change points and athreshold ω_(β) of a change point detection;

at step 2.2, calculating a slope of a degradation segment from aprevious change point to the i-th point; if i−i_(γ)>m, calculating theslope η_(i) of the current segment based on the formula (14), andletting i=i+1:

$\begin{matrix}{{\eta_{i} = \frac{\sum\limits_{j = i_{\gamma}}^{i - 1}\left( {x_{j + 1} - x_{j}} \right)}{i - i_{\gamma} - 1}};} & (14)\end{matrix}$

at step 2.3, calculating a change point detection index β(i) of the i-thpoint based on the formula (15):

$\begin{matrix}{{{\beta (i)} = \frac{\left| {\sum\limits_{j = 1}^{m}\left( {x_{i} + {\eta_{i}j\; \tau} - x_{j}} \right)} \right|}{\sqrt{1 + \eta_{i}^{2}}}};} & (15)\end{matrix}$

at step 2.4, determining whether the change point detection index β(i)of the i-th point exceeds the threshold ω_(β); if β(i)>ω_(β), x_(i)being a change point, and letting i_(γ)=i_(γ)+1 and γ_(i) _(γ) =i; and

at step 2.5, if i≤k, letting i=i+1, and then, performing step 2.2.

Step 11 specifically includes the following steps (with a flow as shownin FIG. 3):

at step 11.1, selecting the number n of Monte Carlo samples, andinitializing parameters i=1 and v_(i,j)=φ_(k);

at step 11.2, generating n random numbers r_(j) obeying a uniformdistribution on [0,1];

at step 11.3, for the j-th Monte Carlo sample sequence, letting i=i+1and v_(i+1,j)=s, where s satisfies

${{\sum\limits_{w = 1}^{s - 1}{p_{v_{i,j},w}(\tau)}} < r_{j} \leq {\sum\limits_{w = 1}^{s}{p_{v_{i,j},w}(\tau)}}},$

and p_(v) _(i,j) _(,w)(τ) refers to a probability that the degradationmode transforms from the v_(i,j)-th mode into the w-th mode over time τ;if i<l_(k)/τ, returning to step 11.2; otherwise, performing step 11.4;

at step 11.4, calculating a total time length that each Monte Carlosample sequence stays in each mode within a time interval (t_(k),t_(k)+l_(k)), and denoting the length as

${S_{s,j} = {\sum\limits_{i = 1}^{l_{k}/\tau}{I\left( {v_{i,j} = s} \right)}}};$

and

at step 11.5, calculating a numerical distribution of ƒ_(Ω(Φ(t) _(k)_(),l) _(k) ₎(x) of Ω[Φ(t_(k)),l_(k)] based on the formula (16):

$\begin{matrix}{{{f_{\Omega {({{\Phi {(t_{k})}},l_{k}})}}(x)} = \frac{\underset{j = 1}{\sum\limits^{n}}{I\left( {{\underset{s = 1}{\sum\limits^{n}}{S_{s,j}{\lambda (s)}}} = x} \right)}}{n}};{when}} & (16) \\{{\sum\limits_{s = 1}^{N}{S_{s,j}{\lambda (s)}}} = x} & \;\end{matrix}$

is established,

${{I\left( {{\sum\limits_{s = 1}^{N}{S_{s,j}{\lambda (s)}}} = x} \right)} = 1};$

otherwise,

${{I\left( {{\sum\limits_{s = 1}^{N}{S_{s,j}{\lambda (s)}}} = x} \right)} = 0}.$

The present disclosure will be described with data of No. 2 blastfurnace in Liuzhou Iron & Steel Company based on a MATLAB tool in thisexample so as to show effects of the present disclosure in combinationwith accompanying drawings.

(1) Temperature sensor data of a furnace wall of a blast furnace iscollected, temperature data (selecting a daily average temperaturevalue) collected continuously by a thermocouple at a height of 8.2meters for 1506 days is selected denoted as x₀, x₁, x₂, . . . , x₁₅₀₆,and degradation data is as shown in FIG. 4.

(2) Change points, denoted as γ₁, γ₂, . . . , before a moment t_(k) aredetected based on a change point detection algorithm, and a change pointdetection result is as shown in FIG. 5 (for example, t_(k)=1300).

(3) A slope of a line segment segmented by the change points iscalculated based on the formula (14) respectively.

(4) A local density ρ_(j) of each line segment and a minimum distanceδ_(j) of a line segment greater than the density of the line segment arecalculated based on the formulas (1), (2) and (3) respectively.

(5) The line segments in which ρ_(j)>2 and δ_(j)>26 are selected asclustering centers, and a total of three clustering centers areobtained, that is, N=3.

(6) The clustering of the line segments is completed by calculating aclustering center closest to each line segment respectively, and aclustering result is obtained as shown in FIG. 6.

(7) An estimation value of a state transition probability matrix ofMarkov chain is obtained based on the formulas (5) and (6).

(8) A estimation value H=0.8959 of a Hurst exponent of a degradationprocess is obtained based on the formula (7).

(9) An estimation value of a drift term coefficient of each segment ofdegradation path is obtained based on the formula (8).

(10) An expectation and a variance of the drift term coefficient in eachdegradation mode are obtained based on the formulas (9) and (10)respectively.

(11) An estimation value σ_(H)=1.0599 of a diffusion term coefficient isobtained based on the formula (11).

(12) A numerical solution ƒ_(Ω[Φ(t) _(k) _(),l) _(k) _(]) of adistribution Ω[Φ(t_(k)),l_(k)] is obtained based on the Monte Carloalgorithm.

(13) If the threshold is 530° C., a prediction result of a remaininguseful life distribution at the k-th sampling moment is obtained basedon the formula (12).

As can be seen from the temperature degradation curve of the furnacewall of the blast furnace wall that the temperature threshold 530° C. ishit for the first time at the 1456-th day. Results shown in FIG. 7 areobtained by performing prediction of the remaining useful life on the1300-th, 1310-th, 1320-th, 1330-th, 1340-th, 1350-th, 1360-th, 1370-th,1380-th, 1390-th and 1400-th days respectively. It can be seen that truevalues of the remaining useful life are all located at positions withhigher probability densities of prediction results, thereby verifyingeffectiveness of the method provided by the present disclosure.

Of course, the above descriptions are not intended to limit the presentdisclosure, and the present disclosure is also not limited to the aboveexamples. Variations, modifications, additions or substitutions made bypersons of ordinary skill in the art shall also belong to the scope ofprotection of the present disclosure.

1. A method of modeling a multi-mode degradation process and predictinga remaining useful life, comprising the following steps: at step 1,collecting degradation data x₀, x₁, x₂, . . . , x_(k) of equipment atequal-interval sampling moments t₀, t₁, t₂, . . . , t_(k) respectively,where a sampling interval is τ, and the number of sampling is k; at step2, detecting slope change points, denoted as γ₁, γ₂, . . . , γ_(j),γ_(j+1) . . . , of a historical degradation process according to achange point detection method; at step 3, obtaining a degradationsegment by taking the points γ_(j) and γ_(j+1) obtained at step 2 asendpoints, calculating a slope η_(γ) _(j+1) , of the degradation segmentbased on the following formula, and taking the slope η_(γ) _(j+1) as afeature value of the j-th degradation segment;$\eta_{\gamma_{j + 1}} = \frac{\underset{j = i_{\gamma}}{\sum\limits^{i_{\gamma + 1}}}\left( {x_{j + 1} - x_{j}} \right)}{i_{\gamma + 1} - 1}$calculating a local density ρ_(j) of the feature value of eachdegradation segment, and calculating a minimum distance δ_(j) of afeature value greater than the local density, wherein the local densityρ_(j) is calculated based on the formula (1): $\begin{matrix}{{\rho_{j} = {\sum\limits_{i}{\chi \left( {d_{ji} - d_{c}} \right)}}},} & (1)\end{matrix}$ wherein d_(c) is a truncation distance,d_(ji)=|η_(i)−η_(j)|, and a function χ(·) is defined as follows:$\begin{matrix}{{\chi (a)} = \left\{ {\begin{matrix}{l,} & {a < 0} \\{0,} & {a \geq 0}\end{matrix};} \right.} & (2)\end{matrix}$ and calculating the minimum distance δ_(j) based on thefollowing formula (3):δ_(j)=min_(i:ρ) _(i) _(>ρ) _(j) d _(ji)  (3); at step 4, if ρ_(j) andδ_(j) are greater than corresponding thresholds respectively, a linesegment obtained with the slope change points γ_(j) and γ_(j+1) asendpoints being one clustering center; according to this method,denoting the number of the obtained clustering centers as N, that is,clustering the line segments segmented by the slope change points as Ncategories according to the slopes of the line segments, denoting adegradation mode of a sampling moment u as Φ(u), letting φ_(i)=Φ(t_(i)),and denoting time points of changes of the degradation mode as c₁, c₂, .. . ; at step 5, establishing a degradation model based on the formula(4):X(t)=X(0)+∫₀ ^(t)λ[Φ(u)]du+σ _(H) B _(H)(t)  (4), where X(0) refers toan initial value of a degradation process, λ[Φ(u)] refers to a driftterm coefficient; when λ(i)˜N(μ_(λ) _(i) , σ_(λ) _(i) ²), σ_(H) refersto a diffusion term coefficient, B_(H) (t) refers to a standard fractalBrownian motion, and the degradation mode Φ(u) is one continuous-timeMarkov chain with a transition probability matrix being Q;${Q = \begin{bmatrix}{- q_{1}} & q_{12} & \ldots & q_{1N} \\q_{21} & {- q_{2}} & \ldots & q_{2N} \\\vdots & \vdots & \ddots & \vdots \\q_{N1} & q_{N2} & \ldots & {- q_{N}}\end{bmatrix}};$ at step 6, estimating q_(j) and q_(ij) in thetransition probability matrix Q of the continuous-time Markov chainaccording to φ₀, φ₁, φ₂, . . . , φ_(k) based on the formulas (5) and(6): $\begin{matrix}{{q_{j} = \frac{m_{j}}{\sum\limits_{i = 1}^{m_{j}}v_{i}^{(j)}}};} & (5) \\{{q_{ij} = {\frac{m_{ij}}{m_{i}}q_{i}}},} & (6)\end{matrix}$ wherein m_(j) refers to a number of times that thedegradation mode hits and stays in the j-th mode before the momentt_(k), m_(ij) refers to a number of times that the degradation modetransits from the i-th mode to the j-th mode before the moment t_(k),and v_(i) ^((j)) refers to a stay time of the degradation mode hittingthe j-th mode at the i-th time; at step 7, estimating a Hurst exponent Hof the degradation process based on the formula (7): $\begin{matrix}{{H = {\frac{1}{2}\log_{2}\frac{E\left\lbrack \left( {\sum\limits_{j = 1}^{p}{\theta_{j}x_{2{({i + j})}}}} \right)^{2} \right\rbrack}{E\left\lbrack \left( {\sum\limits_{j = 1}^{p}{\theta_{j}x_{i + j}}} \right)^{2} \right\rbrack}}},} & (7)\end{matrix}$ wherein θ₁, θ₂, . . . , θ_(p) refer to waveletdecomposition high-pass filter coefficients based on Symlets waveletfunction, p refers to a number of order of a vanishing moment of thewavelet function, and E(·) refers to a mathematic expectation; at step8, estimating an estimation value λ_(i) of a drift term coefficientλ[Φ(u)] in the degradation mode of the i-th segment based on the formula(8) respectively: $\begin{matrix}{{\lambda_{i} = \frac{{\overset{\sim}{x}}_{c_{i}:c_{i + 1}}^{T}Q_{{\overset{\sim}{x}}_{c_{i}:c_{i + 1}}}^{- 1}I_{i}}{\tau \; I_{i}^{T}Q_{{\overset{\sim}{x}}_{c_{i}:c_{i + 1}}}^{- 1}I_{i}}},} & (8)\end{matrix}$ wherein I_(i) refers to one c_(i+1)−c_(i)-dimension columnvector, each element of the column vector is 1, {tilde over (x)}_(c)_(i) _(:c) _(i+1) =[x_(c) _(i) ₊₁−x_(c) _(i) , x_(c) _(i) ₊₂−x_(c) _(i)₊₁, . . . , x_(c) _(i+1) −x_(c) _(i+1) ₊₁]^(T),$Q_{{\overset{\sim}{x}}_{c_{i}:c_{i + 1}}}$ refers to onec_(i+1)−c_(i)-dimension covariance matrix, and the element in the i-throw and the j-th column of the covariance matrix is½[|i−j+1|^(2H)τ^(2H)+|i−j−1|^(2H)τ^(2H)−2|i−j|^(2H)τ^(2H)]; at step 9,estimating an expectation μ_(λ) _(j) and a variance σ_(λ) _(j) of thedrift term coefficient λ[Φ(u)] in each degradation mode based on theformulas (9) and (10) respectively: $\begin{matrix}{{\mu_{\lambda_{j}} = \frac{\sum\limits_{i = 1}^{m_{j}}\lambda_{j,i}}{m_{j}}};} & (9) \\{{\sigma_{\lambda_{j}} = \sqrt{\frac{\sum\limits_{i = 1}^{m_{j}}\left( {\lambda_{j,i} - \mu_{\lambda_{j}}} \right)^{2}}{m_{j}}}},} & (10)\end{matrix}$ wherein λ_(j,i) refers to an estimation value of the driftterm coefficient obtained when the j-th degradation mode is hit at thei-th time; at step 10, estimating a diffusion term coefficient σ_(H)based on the formula (11): $\begin{matrix}{{\sigma_{H} = \sqrt{\frac{\left( {{\overset{\sim}{x}}_{0:k}^{T} - {\Lambda \tau}} \right)^{T}{Q_{{\overset{\sim}{x}}_{0:k}}^{- 1}\left( {{\overset{\sim}{x}}_{0:k}^{T} - {\Lambda\tau}} \right)}}{k}}},} & (11)\end{matrix}$ wherein Λ^(T)=[λ₁I₁ ^(T), λ₂I₂ ^(T), . . . , λ_(m)I_(m)^(T)], {tilde over (x)}_(0:k)=[x₁−x₀, x₂−x₁, . . . , x_(k)−x_(k-1)]^(T),Q_({tilde over (x)}) _(0:k) refers to one k-dimension covariance matrix,and the element in the i-th row and the j-th column of the covariancematrix is ½[|i−j+1|^(2H)τ^(2H)+|i−j−1|^(2H)τ^(2H)−2|i−j|^(2H)τ^(2H)]; atstep 11, letting Ω[Φ(t_(k)),l_(k)]=∫_(t) _(k) ^(t) ^(k) ^(+l) ^(k)λ[Φ(u)]du, and obtaining a numerical distribution ƒ_(Ω[Φ(t) _(k) _(),l)_(k) _(]) of Ω[Φ(t_(k)),l_(k)] by Monte Carlo method; and at step 12,for a given failure threshold ω, an approximate distribution of a firsthitting time of the degradation process being: $\begin{matrix}{{{f_{l}\left( l_{k} \right)} = {\sum\limits_{i = 1}^{N}{\int_{\lambda_{\min}l_{k}}^{\lambda_{\max}l_{k}}{{p_{{\Phi {(t_{k})}}i}\left( l_{k} \right)}{f_{\Omega {({{\Phi {(t_{k})}},t_{k}})}}(s)}\frac{\sigma \left( l_{k} \right)}{\sqrt{2\pi {\sigma (0)}{\int_{0}^{l_{k}}{{\sigma (t)}dt}}}}\left\{ {\frac{\omega - x_{k} - s}{\int_{0}^{l_{k}}{{\sigma (t)}dt}} + \frac{\lambda (i)}{\sigma \left( l_{k} \right)}} \right\} \times \exp \left\{ {- \frac{\left( {\omega - x_{k} - s} \right)^{2}}{2{\sigma (0)}{\int_{0}^{l_{k}}{{\sigma (t)}dt}}}} \right\} {ds}}}}},} & (12)\end{matrix}$ wherein λ_(min)=min{λ(1), λ(2), . . . , λ(N)},λ_(max)=max{λ(1), λ(2), . . . , λ(N)}, p_(Φ(t) _(k)_()i)(l_(k))=P{{|Φ(t_(k)+l_(k))}=i|Φ(t_(k))}, and $\begin{matrix}{{{\sigma (t)} = {\sigma_{H}\begin{Bmatrix}{{\sum\limits_{i = 1}^{\lfloor\frac{t}{\tau}\rfloor}\left\lbrack {\int_{{({i - 1})} \cdot \tau}^{i \cdot \tau}{c_{H}s^{\frac{1}{2}}{\int_{{\lfloor\frac{t}{\tau}\rfloor} \cdot \tau}^{{\lfloor\frac{t + \tau}{\tau}\rfloor} \cdot \tau}{\left( {u - s} \right)^{H - \frac{3}{2}}u^{H - \frac{1}{2}}{duds}}}}} \right\rbrack^{2}} +} \\\left\lbrack {\int_{{\lfloor\frac{t}{\tau}\rfloor} \cdot \tau}^{{\lfloor\frac{t + \tau}{\tau}\rfloor} \cdot \tau}{c_{H}s^{\frac{1}{2}}{\int_{s}^{{\lfloor\frac{t + \tau}{\tau}\rfloor} \cdot \tau}{\left( {u - s} \right)^{H - \frac{3}{2}}u^{H - \frac{1}{2}}{duds}}}}} \right\rbrack^{2}\end{Bmatrix}^{\frac{1}{2}}}};} & \left. 13 \right) \\{\mspace{79mu} {{c_{H} = \sqrt{\frac{2H{\Gamma \left( {\frac{3}{2} - H} \right)}}{{\Gamma \left( {\frac{1}{2} + H} \right)}{\Gamma \left( {2 - {2H}} \right)}}}},}} & \;\end{matrix}$ wherein Γ(·) refers to a Gamma function.
 2. The methodaccording to claim 1, wherein the step 2 specifically comprises thefollowing steps: at step 2.1, initializing parameters, letting γ₁=0,i=1, and i_(γ)=1, and selecting a minimum interval mτ between two changepoints and a threshold ω_(β) of a change point detection; at step 2.2,calculating the slope of the degradation segment from a previous changepoint to the i-th point; when i−i_(γ)>m, calculating the slope η_(i) ofthe current segment based on the formula (14), and letting i=i+1:$\begin{matrix}{{\eta_{i} = \frac{\sum\limits_{j = i_{\gamma}}^{\;_{i - 1}}\left( {x_{j + 1} - x_{j}} \right)}{i - i_{\gamma} - 1}};} & (14)\end{matrix}$ at step 2.3, calculating a change point detection indexβ(i) of the i-th point based on the formula (15): $\begin{matrix}{{{\beta (i)} = \frac{{\sum\limits_{j = 1}^{m}\left( {x_{i} + {\eta_{i}j\; \tau} - x_{j}} \right)}}{\sqrt{1 + \eta_{i}^{2}}}};} & (15)\end{matrix}$ at step 2.4, determining whether the change pointdetection index β(i) of the i-th point exceeds the threshold ω_(β); whenβ(i)>ω_(β), x_(i) being a change point, and letting i_(γ)=i_(γ)+1 andγ_(i) _(γ) =i; and at step 2.5, when i≤k, letting i=i+1, and thenperforming step 2.2.
 3. The method according to claim 1, wherein thestep 11 specifically comprises the following steps: at step 11.1,selecting the number n of Monte Carlo samples, and initializingparameters i=1 and v_(i,j)=φ_(k); at step 11.2, generating n randomnumbers r_(j) obeying a uniform distribution on [0,1]; and at step 11.3,for the j-th Monte Carlo sample sequence, letting i=i+1 and v_(i+1,j)=s,wherein s satisfies${{\sum\limits_{w = 1}^{s - 1}{p_{v_{i,j},w}(\tau)}} < r_{j} \leq {\sum\limits_{w = 1}^{s}{p_{v_{i,j},w}(\tau)}}},$and p_(v) _(i,j) _(,w)(τ) refers to a probability that the degradationmode transforms from the v_(i,j)-th mode into the w-th mode over time τ;when i<l_(k)/τ, returning to step 11.2; otherwise, performing step 11.4;at step 11.4, calculating a total time length of each Monte Carlo samplesequence staying in each mode within a time interval (t_(k),t_(k)+l_(k)), and denoting the length as${S_{s,j} = {\sum\limits_{i = 1}^{l_{k}/\tau}{I\left( {v_{i,j} = s} \right)}}};$and at step 11.5, calculating a numerical distribution ƒ_(Ω(Φ(t) _(k)_(),l) _(k) ₎(x) of Ω[Φ(t_(k)),l_(k)] based on the formula (16):$\begin{matrix}{{{f_{\Omega {({{\Phi {(t_{k})}},l_{k}})}}(x)} = \frac{\sum\limits_{j = 1}^{n}{I\left( {{\sum\limits_{s = 1}^{N}{S_{s,j}{\lambda (s)}}} = x} \right)}}{n}};} & (16) \\{{{when}\mspace{14mu} {\sum\limits_{s = 1}^{N}{S_{s,j}{\lambda (s)}}}} = x} & \;\end{matrix}$ is established,${{I\left( {{\sum\limits_{s = 1}^{N}{S_{s,j}{\lambda (s)}}} = x} \right)} = 1};$otherwise,${I\left( {{\sum\limits_{s = 1}^{N}{S_{s,j}{\lambda (s)}}} = x} \right)} = 0.$